######################### Simulations
pacman::p_load("riclpmr","dagitty")
set.seed(1234)
RICLPM5 <- riclpm_text(var_groups = list(x=c("x1","x2","x3"),y=c("y1","y2","y3")))
CLPM_u <- '
  # Estimate the lagged effects between the observed variables.
  x2 + y2 ~ x1 + y1
  x3 + y3 ~ x2 + y2
  x3 ~~ x3
  x2 ~~ x2
  x1 ~~ x1
  y3 ~~ y3
  y2 ~~ y2
  y1 ~~ y1
  y3 ~~ x3
  y2 ~~ x2
  y1 ~~ x1
'
## no contemporaneous effects
g <- dagitty('dag{
                  y1 -> x1
                  x1 -> y1
                  x1 -> x2 [beta=.7]
                  y1 -> y2 [beta=.7]
                  x1 -> y2 [beta=.2]
                  y1 -> x2 [beta=.2]
                  y2 -> x2 
                  x2 -> y2
                  y2 -> y3 [beta=.7]
                  x2 -> y3 [beta=.2]
                  x2 -> x3 [beta=.7]
                  y2 -> x3 [beta=.2]
                  y3 -> x3
                  x3 -> y3
 }')

plot(g)
library(retry)


riclpm_uncorr <- list()
for(i in 1:1000){
  x <- retry(simulateSEM( g,b.default = 0,N = 1000),when = "'Sigma' is not positive definite")
  sim_ri <-
    riclpmr::lavriclpm(
      RICLPM5,
      data = x)
  riclpm_uncorr[[i]] <- parameterestimates(sim_ri) %>% filter(label %in% c("x_y","y_x")) %>% slice(3)
}

clpm_uncorr <- list()
for(i in 1:1000){
  x <- retry(simulateSEM( g,b.lower = .1,b.upper=.3,N = 1000),when = "'Sigma' is not positive definite")
  sim_ri <-
    sem(
      CLPM_u,
      data = x)
  clpm_uncorr[[i]] <- parameterestimates(sim_ri) %>% slice(3)
}


g <- dagitty('dag{
                  y1 -> x1
                  x1 -> y1
                  x1 -> x2 [beta=.7]
                  y1 -> y2 [beta=.7]
                  x1 -> y2 [beta=.2]
                  y1 -> x2 [beta=.2]
                  y2 -> x2 
                  x2 -> y2
                  y2 -> y3 [beta=.7]
                  x2 -> y3 [beta=.2]
                  x2 -> x3 [beta=.7]
                  y2 -> x3 [beta=.2]
                  y3 -> x3
                  x3 -> y3
 }')
plot(g)

riclpm_corr <- list()
for(i in 1:1000){
  x <- retry(simulateSEM( g,b.lower = .1,b.upper=.3,N = 1000),when = "'Sigma' is not positive definite")
  sim_ri <-
    riclpmr::lavriclpm(
      RICLPM5,
      data = x)
  riclpm_corr[[i]] <- parameterestimates(sim_ri) %>% filter(label %in% c("y_x")) %>% slice(1)
}
library(lavaan)

clpm_corr <- list()

for(i in 1:1000){
  x <- retry(simulateSEM( g,b.lower = .1,b.upper=.3,N = 1000),when = "'Sigma' is not positive definite")
  sim_ri <-
    sem(
      CLPM_u,
      data = x)
  clpm_corr[[i]] <- parameterestimates(sim_ri) %>% slice(2,3,6,7)
}

simresults <-
  rbind(
    data.frame(est=data.table::rbindlist(riclpm_uncorr)$est,model="RI-CLPM",contemporaneous="no"),
    data.frame(est=data.table::rbindlist(clpm_uncorr)$est,model="CLPM",contemporaneous="no"),
    data.frame(est=data.table::rbindlist(clpm_corr)$est,model="CLPM",contemporaneous="yes"),
    data.frame(est=data.table::rbindlist(riclpm_corr)$est,model="RI-CLPM",contemporaneous="yes"))

simresults %>% group_by(model,contemporaneous) %>% summarise(mean(est))
p1 <- ggplot(simresults,aes(x=model,fill=contemporaneous,y=est))+geom_boxplot() + geom_hline(yintercept = .2) + ylim(-.5,.5) + theme_bw()
ggsave("tables/simulation.pdf",width=10,height=8)

(0-.20)/.20
(.30-.20)/.20
